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The  Gibbs  free  energy  of  solvation  and  dissociation  of  hydrogen  chloride  in  water  is  calculated  through 
a  combined  molecular  simulation/quantum  chemical  approach  at  four  temperatures  between  T  =  300 
and  450  K.  The  Gibbs  free  energy  is  first  decomposed  into  the  sum  of  two  components:  the  Gibbs  free 
energy  of  transfer  of  molecular  HCI  from  the  vapor  to  the  aqueous  liquid  phase  and  the  standard-state 
Gibbs  free  energy  of  acid  dissociation  of  HCI  in  aqueous  solution.  The  former  quantity  is  calculated 
using  Gibbs  ensemble  Monte  Carlo  simulations  using  either  Kohn-Sham  density  functional  theory  or  a 
molecular  mechanics  force  field  to  determine  the  system's  potential  energy.  The  latter  Gibbs  free  energy 
contribution  is  computed  using  a  continuum  solvation  model  utilizing  either  experimental  reference 
data  or  micro-solvated  clusters.  The  predicted  combined  solvation  and  dissociation  Gibbs  free  energies 
agree  very  well  with  available  experimental  data. 


1  Introduction 

The  interaction  between  hydrogen  chloride  and  water  has  been 
studied  for  many  years  due  to  the  importance  of  HCI  as  a  strong 
acid  (including  acid  dissociation  in  very  cold  nanoclusters1) 
and  of  HCl-ice  systems  and  chloride  ions  at  the  air/water 
interface  in  atmospheric  chemistry.2-5  Interest  has  also  been 
shown  in  the  supercritical  water-HCl  system,  as  supercritical 
water  is  becoming  increasingly  popular  for  applications  such  as 
waste  disposal;  the  nature  of  supercritical  water,  however, 
drastically  alters  the  behavior  of  strong  acids.6 
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Studies  of  the  thermophysical  properties  of  HCI  dissolution 
in  water  were  performed  early  in  the  last  century  using  electro¬ 
chemical  cells,7’8  and  the  structures  of  both  dilute  and  concen¬ 
trated  HCI  solutions  have  been  probed  in  neutron  and  X-ray 
diffraction  experiments,9’10  including  studies  on  the  chloride  ion 
solvation  structure.11  Experimental  measurements  of  the  vapor 
pressures  of  HCI  over  water  have  been  reported  for  various 
concentrations  and  temperatures,12-19  including  vapor-  and 
liquid-phase  mole  fractions  of  HCI  as  functions  of  tempera¬ 
ture.20  Equations  to  predict  the  HCI  vapor  pressure,21-25  as  well 
as  its  Henry’s  law  constants,26  have  also  been  developed  based 
on  previous  experimental  data. 

There  is  some  disagreement  in  the  literature  concerning  the 
value  of  the  acid  dissociation  constant  of  HCI  in  aqueous 
solution  {K^ aq).  Ruaya  and  Seward27  report  values  for  a  wide 
range  of  temperatures,  including  a  value  of  0.71  at  T  =  298  K, 
showing  good  agreement  with  some  previous  studies.  However, 
a  value  of  /Cv,q  less  than  unity  would  imply  a  significant 
amount  of  undissociated  HCI  at  room  temperature,  which  is 
unexpected  at  low  temperatures  for  a  strong  acid,  and  Ruaya 
and  Seward27  mention  other  values  that  are  about  an  order  of 
magnitude  larger.  Much  larger  values  make  more  intuitive  sense 
and  have  been  reported  by  various  research  groups.23’24’28’29  It 
is  also  of  interest  to  note  that  Clegg  and  Brimblecombe24 
recommended  using  a  form  of  Hemy’s  law  constant  that  avoids 
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the  need  to  know  “inaccessible  parameters”  including  the  acid 
dissociation  constant,  and  commented  on  the  difficulty  in  deter¬ 
mining  its  true  value.  Balbuena  et  al .30  report  the  Gibbs  free  energy 
of  HC1  dissociation  in  water  along  the  water  saturation  curve  below 
the  critical  temperature,  ranging  from  about  —11  kcal  mol-1  at 
T  =  300  K  to  around  —5  kcal  mol  1  near  T  =  500  I<  (based  on 
experimentally  measured  acid  dissociation  constants31).  Ho  et  al.32 
give  a  formula  to  estimate  at  elevated  temperatures  as  a 
function  of  the  solvent  density  and  temperature.  At  T  =  473  I< 
(the  lowest  temperature  at  which  their  formula  is  valid),  this 
formula  predicts  log (Ka>aq)  =  0.55,  which  gives  a  Gibbs  free  energy 
of  acid  dissociation,  AGaaq,  of  1.2  kcal  mol-1.  Robinson,28 
Marsh  and  McElroy,23  and  Pokrovskii33  all  give  values  of  AGa  aq  = 
—8.4  kcal  mol-1  at  T  =  298  K.  Pokrovskii33  points  out  a  discrepancy 
between  their  values  and  those  of  Ruaya  and  Seward  27  and  argues 
that  the  data  used  by  the  latter  group  can  be  interpreted  in 
multiple  ways. 

In  contrast,  the  Henry’s  law  constant  of  HCl  in  water  appears  to 
be  less  ambiguous.  Clegg  and  Brimblecombe34  report  a  value  of 
Kh  =  1.84  x  106  mol2  kg-2  atm-1  under  ambient  conditions.  The 
following  year,  Clegg  and  Brimblecombe24  reported  values  around 
2.0  x  106  mol2  kg-2  atm-1  for  concentrations  of  HCl  in  solution 
between  0  and  200  mol  kg-1  by  using  data  from  a  variety  of 
sources,  while  Brimblecombe  and  Clegg26  recommend  a  value  of 
2.04  x  106  mol2  kg-2  atm-1  at  T  =  298  K.  Marsh  and  McElroy23 
calculated  equilibrium  constants  for  both  the  transfer  of  mole¬ 
cular  HCl  from  the  vapor  to  aqueous  phase  (dimensionless)  and 
the  dissociation  of  HCl  in  water  (in  units  of  mole  dm-3).  Combin¬ 
ing  these  values  ( via  eqn  (11)  in  their  paper)  and  converting  units 
results  in  /CH  =  1.85  x  106  mol2  kg-2  atm-1,  which  is  very  close  to 
previous  values,  although  Clegg  and  Brimblecombe24  note  that 
the  disagreement  grows  at  lower  temperatures. 

The  simulation  literature  surrounding  the  HCl-water  system  is 
just  as  extensive;  in  the  interest  of  brevity,  the  following  summary 
will  only  be  concerned  with  calculations  involving  the  bulk  aqueous 
system  and  bulk  solvation  {i.e.,  studies  dealing  only  with  clusters  or 
ice  will  be  ignored).  Molecular  simulations  of  hydrogen  chloride  in 
water  have  mostly  been  performed  with  the  molecular  dynamics 
technique,  including  simulations  from  first-principles35-41  and 
using  empirical  force  fields, 30,42-45  including  a  study  of  HCl  at 
the  air-water  interface  using  QM/MM  techniques.46  Some  of  these 
studies  compute  the  p/Qjilq  of  HCl  in  water,30’42’43  others  compute 
the  potential  of  mean  force,44  while  the  rest  examine  the  structure 
of  the  solvating  water  molecules  around  the  ions.  A  few  of  the 
simulations  were  performed  at  elevated  temperatures,30’37’44  mostly 
near  the  critical  temperature  of  water. 

Less  frequently,  Monte  Carlo  approaches  have  also  been  used  to 
examine  the  aqueous  HCl  system,  although  until  this  point  all 
calculations  employing  this  method  have  used  empirical  force 
fields  to  represent  the  solvated  HCl  molecule;  this  includes  some 
work  on  the  Gibbs  free  energy  curves  along  the  reaction  coordinate 
to  determine  the  most  probable  dissociation  mechanism  47  and 
other  work  using  Monte  Carlo  simulations  to  reinterpret  previous 
experimental  data  to  determine  the  structure  48,49 

In  addition  to  molecular  simulation,  the  Gibbs  free  energy  of 
solvation  of  hydrogen  chloride  in  water  has  been  computed  with 


various  other  techniques.  Among  these  are  QM/MM  methods50’51 
and  calculation  from  experimental  data,52  such  as  gas-phase 
proton  affinities  and  aqueous  p KaAq  values.53  The  structure  of 
concentrated  HCl  solutions  has  also  been  determined  by  empiri¬ 
cally  fitting  X-ray  scattering  data.54 

In  the  above  molecular  simulation  studies,  no  effort  was 
made  to  compute  the  Gibbs  free  energy  of  aqueous  solvation  for 
HCl  or  its  Henry’s  law  constant.  The  goal  of  the  current  study  is 
to  compute  these  thermodynamic  properties  across  a  tempera¬ 
ture  range  of  150  K  (T  =  300-450  K),  through  the  use  of  first 
principles  Monte  Carlo  simulations  and  quantum  mechanical 
continuum  solvation  approaches.  The  chemical  equilibria  for 
the  solvation  and  dissociation  of  gaseous  HCl  are  given  by: 


HCl(g)  <->  H(aq)+  +  Cl(aq) 

(1) 

HCl(g)  <->  HCl(aq) 

(2) 

HCl(aq)  <->  H(aq)  1  +  Cl(aq) 

(3) 

law  constant  of 

„  cm  ca-  2 

Kh  —  7  hci 

(4) 

Phc\ 


where  cx  is  the  concentration  of  species  X  in  the  aqueous  phase, 
yHci  is  the  geometric  mean  activity  coefficient  of  H+and  Cl-  in 
solution,  and  /jh,  i  is  the  partial  pressure  of  HCl  over  the 
solution.26  It  should  be  noted  that  other  definitions  of  KH  can 
be  found  in  the  literature.  Eqn  (4)  can  be  rewritten  in  terms  of 
the  Gibbs  free  energy  changes  associated  with  eqn  (2)  and  (3): 


Aji  =  exp 


AGcomb  \  (c9)2 
RT  )  p° 


(5) 


AGcomb  AGtrans  +  A  Gaaq  (6) 

where  AGttans,  AG°aq,  and  AGcomb  are  the  Gibbs  free  energy  of 
transfer  of  molecular  HCl  (obtained  from  the  ratio  of  number 
densities),  the  standard-state  Gibbs  free  energy  of  acid  dissociation 
in  aqueous  solution,  and  the  Gibbs  free  energy  of  the  combined 
solvation-acid  dissociation  process,  respectively.  and  ce  are  the 
standard  pressure  and  the  standard  concentration,  respectively. 

In  this  work,  particle-based  Monte  Carlo  simulations  are 
carried  out  to  compute  AGtrans,  and  continuum  solvation  models 
are  used  to  compute  AGa  aq.  Notice  that  AGtrans  is  not  equivalent  to 
the  standard-state  Gibbs  free  energy  of  transfer  as  the  simulations 
are  carried  out  at  the  equilibrium  vapor  pressure  of  HCl  above  the 
aqueous  solution  (in  the  case  of  the  first-principles  simulations) 
or  at  the  equilibrium  vapor  pressure  of  the  binary  mixture  (for  the 
simulations  with  an  empirical  force  field).  Further  details  are 
given  below,  but  in  both  cases  the  pressure  is  a  function  of 
the  temperature.  The  standard-state  Gibbs  free  energy  of  acid 
dissociation  of  HCl  in  aqueous  solution  is  given  by: 

AG°aq  =  GV+Gaq,a—  G°q,HCl  (?) 

and  Gaq>x  is  defined  as 

Gaq,x  =  Gg,x  +  AGsiX  (8) 
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The  first  term  in  eqn  (8)  is  the  standard-state  Gibbs  free  energy 
of  the  species  in  an  ideal  gas  at  1  bar  pressure,  and  the  second 
term  is  the  standard-state  Gibbs  free  energy  of  solvation,  i.e., 
the  Gibbs  free  energy  of  transfer  of  a  species  from  the  gas  phase 
at  a  solute  partial  pressure  of  1  bar  to  a  1  M  ideal  solution.  The 
standard-state  Gibbs  free  energy  of  solvation  is  expressed  as 

RT 

A^,x  =  AG^x  +  /mn—  (9) 

where  the  first  term  is  the  fixed-concentration  Gibbs  free  energy  of 
solvation,  in  other  words,  the  Gibbs  free  energy  of  transfer  of  the 
solute  from  the  gas  phase  with  a  solute  concentration  of  1  mol  L  1 
to  1  M  solution.  Here,  we  take  V*  =  1  L  mol  1  and  p  =  1  bar. 

It  is  important  to  emphasize  some  notational  conventions. 
AGs,X)  the  standard-state  Gibbs  free  energy  of  solvation,  and 
AGtransj  the  Gibbs  free  energy  of  transfer,  describe  the  transfer 
process  of  a  molecule-ion  from  the  gas  phase  to  the  liquid 
phase.  However,  they  are  labeled  distinctly  due  to  the  fact  that 
AGs,x  is  always  computed  for  the  case  of  infinitely  dilute  gases 
and  solutions  and  converted  to  standard-state  conditions  with 
the  convention  of  an  ideal  gas  at  P  =  1  bar  and  an  ideal  solution 
at  c  =  1  M,  while  AG^^s  directly  relates  to  an  equilibrium  constant 
and  involves  no  ideality  convention.  In  this  manuscript,  the 
quantities  are  also  computed  with  different  methods,  making  such 
a  distinction  even  more  useful.  One  similar  quantity,  AGcnnlh, 
indicates  the  combined  process  of  the  Gibbs  free  energy  of  transfer 
of  HC1  and  its  dissociation  in  the  solvent  (see  eqn  (6)). 

The  layout  of  the  remainder  of  this  paper  is  as  follows. 
Section  2  describes  the  details  of  the  particle-based  Monte 
Carlo  simulations  and  of  the  different  methods  for  computing 
the  standard-state  Gibbs  free  energy  of  acid  dissociation  with 
an  implicit  solvation  model.  Section  3  presents  the  resulting 
Gibbs  free  energies  of  transfer  from  the  vapor  to  the  liquid 
phase  at  several  temperatures,  combines  these  with  the  standard- 
state  Gibbs  free  energy  of  acid  dissociation  calculations  to 
estimate  the  Gibbs  free  energy  of  the  combined  solvation-acid 
dissociation  of  HCl  in  water,  and  compares  the  results  to 
experimental  results  in  the  literature.  Section  4  summarizes 
the  findings  reported  here. 

2  Computational  details 

2.1  Gibbs  free  energy  of  transfer 

2.1.1  First  principles  Monte  Carlo  simulations.  In  order  to 
calculate  the  Gibbs  free  energy  of  transfer  of  molecular  HCl 
from  the  gas  to  the  aqueous  phase,  first  principles  Monte  Carlo 
(FPMC)  simulations55’56  were  run  in  the  NVT-Gibbs  ensemble57’58 
using  the  CP2K  software.59’60  The  simulated  system  consisted 
of  one  HCl  and  63  H20  molecules  in  two  boxes.  Due  to  the 
expense  of  these  simulations,  water  molecules  were  confined  to 
the  liquid-phase  simulation  box,  i.e.,  the  gas  phase  is  assumed 
to  be  ideal  since  the  single  HCl  molecule  does  not  interact 
with  any  other  molecules.  NVT-Gibbs  ensemble  simulations 
require  particle  transfer  moves  between  the  phases  to  equili¬ 
brate  the  chemical  potential  of  a  given  species,  volume  moves 
to  reach  mechanical  equilibrium,  and  translational,  rotational, 
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and  conformational  moves  to  reach  thermal  equilibrium.  In  the 
present  FPMC  simulations,  the  probabilities  of  performing  a 
given  move  type  were  as  follows:  50%  to  swap  HCl  between  the 
boxes,  5%  to  change  the  volume  of  the  boxes,  and  the  remaining 
45%  was  divided  equally  between  molecular  translations, 
rotations,  and  conformational  changes  of  water  molecules 
(with  the  bond  length  of  HCl  being  held  rigid).  Because  of 
the  larger  fraction  of  water  molecules,  95%  of  the  rotations  and 
translations  were  attempted  on  water.  Maximum  displace¬ 
ments  for  molecular  translations,  rotations  around  the  center 
of  mass,  conformational  changes,  and  volume  moves  were 
adjusted  during  the  equilibration  period  to  give  acceptance 
rates  of  »  50%.  Four  independent  simulations  were  carried  out 
at  each  of  four  temperatures  ( T=  300,  350,  400,  and  450  K),  and 
the  production  periods  consisted  of  1000  Monte  Carlo  cycles, 
where  each  cycle  consists  of  64  FPMC  moves. 

In  these  FPMC  simulations,  the  potential  energy  of  the 
interacting  system  was  computed  using  Kohn-Sham  density 
functional  theory  (KS-DFT)61  with  the  Becke-Lee-Yang-Parr 
(BLYP)  exchange  and  correlation  functionals,62’63  a  triple-^  basis 
set  with  two  polarization  functions,  and  the  norm-conserving 
Goedeclcer-Teter-Hutter  pseudopotentials.64  Use  of  a  reference  cell 
for  the  plane-wave  grids  allowed  the  use  of  a  relatively  low  plane- 
wave  cutoff  of  280  Ry  for  the  electronic  density  that  has  the  benefit 
of  giving  fairly  accurate  liquid  densities  and  Gibbs  free  energies  of 
transfer  for  neat  water.56,65’66  Previous  simulation  studies  have  also 
shown  qualitative  agreement  between  plane-wave  BLYP  simula¬ 
tions  and  higher-level  methods  on  small  H20-HC1  clusters  at 
very  low  temperatures.67’68 

As  in  previous  FPMC  simulations,55’56  empirical  biasing 
potentials  are  used  for  pre-sampling  sequences69’70  including 
molecular  translations,  rotations,  and  conformational  changes 
within  a  single  FPMC  move  and  for  configurational-bias  swap 
moves.71-73  The  empirical  biasing  potentials  were  parameterized 
to  reproduce  the  energy  differences  between  configurations 
obtained  with  KS-DFT.  These  biasing  potentials  used  Lennard- 
Jones  (only  on  O  and  Cl)  and  Coulombic  interactions  of  atomic 
sites  for  the  intermolecular  interactions  [s/kB  =  120  K  for  Cl  and 
89.3  K  for  O,  a  =  3.38  A  for  Cl  and  3.19  A  for  O,  q  =  —0.15  |e|  for  Cl 
and  —0.784  |e|  for  O).  The  Lorentz-Berthelot  combining  rules74 
were  used  for  unlike  interactions,  and  all  interactions  were 
truncated  and  shifted  to  zero  at  rcut  =  6.0  A.  Intramolecular 
parameters  for  harmonic  bond  stretching  and  bending  were  also 
fit  for  water  (rOH  =  0.983  A,  krlkB  =  1.181  x  106  IC  A-2,  0HOH  = 
101.2°,  and  k0lkB  =  15.07  K  deg-2),  while  HCl  was  kept  rigid 
(fHci  =  1-35  A)  throughout  the  simulations.  To  improve  the 
statistics  for  the  computation  of  AGtrans,  a  balancing  factor75’76 
was  added  that  changes  the  potential  energy  of  HCl  in  the  vapor 
phase  and  allows  one  to  find  HCl  with  appreciable  probabilities 
in  both  phases  throughout  a  simulation.  This  balancing  factor 
was  determined  during  a  short  pre-simulation  (A qlkB  =  3600  K 
at  T-  300  K  and  2700  K  for  the  other  three  temperatures),  and  it 
is  removed  in  the  calculation  of  AGtrans  from  the  ratio  of 
number  densities.77 

2.1.2  Monte  Carlo  simulations  with  empirical  force  field. 

Gibbs  ensemble  Monte  Carlo  simulations  were  also  run  with  an 
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empirical  force  field  (EFF)  at  the  same  four  temperatures  to 
assess  the  system-size  dependence  and  to  provide  an  additional 
data  set.  A  two-site  HCl  model  with  Lennard-Jones  and  Coulomb 
potentials  was  fit  to  reproduce  the  experimental  critical  tem¬ 
perature  and  liquid  density  near  the  triple  point  (sa/kB  =  175  K, 
del  ~  3.5  A,  qc  1  =  —0.17  |e|,  eH/kB  =  40  K,  and  aH  =  2.1  A).  The 
TIP4P  model78  was  used  for  water.  The  systems  consisted  of 
one  HCl  molecule  and  either  64,  128,  or  500  water  molecules. 
Approximately  5  x  104  cycles  were  used  for  equilibration,  with 
another  2  x  105  cycles  for  production  (where  a  cycle  is  N  steps, 
with  N  being  the  total  number  of  molecules  in  the  system).  In 
these  simulations,  the  water  molecules  were  also  allowed  to 
transfer  between  the  two  phases. 

2.2  Standard-state  Gibbs  free  energy  of  acid  dissociation 

In  order  to  compute  AGa,aq  via  eqn  (7),  two  different  methods 
were  employed  here.  Method  SM-ref  calculates  A G°  aq  at  different 
temperatures  using  experimental  reference  values  for  the  gas- 
phase  standard-state  Gibbs  free  energy  of  acid  dissociation  and 
for  the  solvation  Gibbs  free  energies  of  H+  and  CP  at  298  K. 
Method  SM-cluster  does  not  require  any  reference  data  but 
employs  micro-solvated  clusters  in  the  continuum  solvation 
model  calculations. 

2.2.1  Method  SM-ref.  One  can  rewrite  eqn  (7)  as 

AG“aq  =  AG°g  +  AG°  H+  +  AGg ,  cr  -  AG°  HC1  (10) 

where  the  first  term  on  the  right  side  is  the  standard-state  Gibbs 
free  energy  of  acid  dissociation  of  HCl  in  the  gas  phase  at  1  bar. 
The  temperature  dependence  of  this  term  is  calculated  by: 

A G°g(T)  =  AG°,g,ref(298)  +  AG°g,caic(T)  -  AG°g,calc(298) 

(11) 

where  the  first  term  is  the  experimental  standard-state  gas- 
phase  Gibbs  free  energy  of  acid  dissociation  at  298  K.79  The 
second  and  third  terms  are  obtained  from  the  Ga  g  of  individual 
compounds  calculated  by  Gaussian0980  at  the  gas/M06-2X/ 
MG3S  level  of  theory.81-83 

The  standard-state  Gibbs  free  energy  of  solvation  of  HCl 
at  different  temperatures  is  obtained  by  eqn  (9)  from  the  fixed- 
concentration  Gibbs  free  energy  of  solvation,  AG|  Ha(r),  calculated 
by  the  SM8T  solvation  model84’85  at  the  SM8T/M06-2X/6-31G(d)//gas/ 
M06-2X/MG3S  level  of  theoiy  using  a  locally  modified  version  of 
Gaussian  09.  Note  that  SM8T  at  298  K  is  identical  to  SM8.86  The 
standard-state  Gibbs  free  energies  of  solvation  of  the  proton  and 
chloride  anion  are  obtained  from  the  corresponding  fixed- 
concentration  Gibbs  free  energies  of  solvation  calculated  by 
the  following  equation 

=  AGgXref(298)  +  AGgXcalc(r)  —  AGgXcalc(298) 

(12) 

where  X  =  H+  and  Cl-.  The  first  term  in  eqn  (12)  is  the  reference 
fixed-concentration  solvation  Gibbs  free  energies  of  the  proton 
(—265.9  kcal  mol-1)  and  the  chloride  anion  (—74.6  kcal  mol-1)  at 
T  =  298  K  obtained  using  data  from  Tissandier  et  al.87  The  last 
two  terms  in  eqn  (12)  are  calculated  by  SM8T/M06-2X/6-31G(d). 
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2.2.2  Method  SM-cluster.  The  standard  state  Gibbs  free 
energy  of  solute  X  in  aqueous  solution,  Gaq  X,  can  be  computed 
as  follows 

^aq,X  =  ^aq,X)iX  ~  n^aq,H20  ~  «X  RT  ln(cwal)  (13) 

where  the  notation  X«x  denotes  a  supersolute  with  nx  water 
molecules  added  explicitly  to  the  molecule  X.  The  last  term  of 
eqn  (13)  is  the  concentration  term,  equal  to  «  2.38  kcal  mol-1  for 
nx  =  1  at  298  K.  The  quantity  cwat  =  pwat/Mvat  is  the  concentration 
of  liquid  water  in  mol  L-1  (cwat  x  55.3  M  at  298  I<). 

We  used  nx  =  0,  1,  and  3  for  HCl,  Cl-,  and  H+,  respectively, 
which  is  equivalent  to  studying  the  reaction 

HCl  +  4H20  o  H703+  +  H2O  Cl-  (14) 

The  quantities  Gaq  on  the  right  hand  side  of  eqn  (13)  were 
calculated  via  eqn  (8),  using  the  GgjX  energies  of  individual 
compounds  computed  with  Gaussian09  by  the  gas/M06-2X/ 
MG3S  method  and  the  fixed-concentration  Gibbs  free  energies 
of  solvation  ^AGgX^  obtained  from  SM8T/M06-2X/6-31G(d)// 
gas/M06-2X/MG3S. 

It  is  difficult  to  estimate  the  reliability  of  eqn  (8)  with  the 
SM8  solvation  model  for  clusters,  but  we  note  that  previous 
tests88’89  of  similar  models  for  Gibbs  free  energies  of  aqueous 
solvation  of  clustered  ions  gave  mean  unsigned  errors  in  the 
range  3. 3-5. 3  kcal  mol-1.  Note  though  that  eqn  (10)  involves 
both  a  proton  and  a  chloride.  In  such  a  case  the  proton  and 
chloride  hydration  Gibbs  free  energies  could  be  combined 
together  to  give  the  hydration  Gibbs  free  energy  of  a  fully 
dissociated  HCl  molecule  (at  infinite  dilution).  This  would 
not  have  any  material  effect  on  the  resulting  total  numbers, 
but  it  would  avoid  references  to  estimated  single-ion  quantities 
and  avoid  the  portion  of  the  error  associated  with  separating 
Gibbs  free  energy  values  pertaining  to  the  neutral  species  into 
contributions  from  oppositely  charged  ions  and  with  the  unknown 
surface  potential  of  water  and  its  temperature  dependence  (the 
surface  contributions  to  cations  and  anions  will  cancel  out  for 
neutral  species,  whether  or  not  they  are  dissociated).  Therefore 
the  errors  might  actually  be  less  than  those  estimated  based  on 
single-ion  values. 

3  Results  and  discussion 

The  Gibbs  free  energies  of  transfer  for  molecular  HCl  deter¬ 
mined  from  the  Gibbs  ensemble  simulations  with  the  BLYP 
functionals  or  the  EFF  are  summarized  in  Table  1.  The  uncer¬ 
tainties  for  the  Monte  Carlo  EFF  and  BLYP  simulations  are 
about  0.1  and  1  kcal  mol-1,  respectively.  These  estimates  are 
the  standard  errors  of  the  mean  taken  from  four  independent 
simulations  for  each  model-temperature  combination.  Values 
for  EFF/64  and  BLYP/63  (where  the  number  after  the  solidus  is 
the  number  of  water  molecules  in  the  simulation)  at  T  =  300  K 
are  not  provided  because  for  the  former  the  simulation  box  for 
the  liquid  phase  became  smaller  than  allowed  by  the  potential 
truncation  and  for  the  latter  the  number  of  successful  swap 
moves  was  deemed  too  small.  Even  at  the  three  higher 
temperatures,  long  periods  without  successful  swap  moves  were 
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Table  1  Gibbs  free  energies  of  transfer  (in  kcal  mol-1)  for  molecular  HCI  from 
the  vapor  phase  to  aqueous  solution  computed  using  the  BLYP  functionals  and 
the  EFF  for  different  system  sizes 


Model/NH  o 

T[  K] 

BLYP/63 

EFF/64 

EFF/128 

EFF/500 

300 

N/A 

N/A 

0.6 

0.8 

350 

1.9 

1.1 

1.0 

1.0 

400 

2.9 

1.2 

1.2 

1.1 

450 

3.6 

1.0 

1.2 

1.2 

observed,  but  these  were  interspersed  by  periods  with  relatively 
high  frequencies  of  particle  transfers.  Analysis  of  the  trajectories 
showed  that  the  HCI  molecule  can  become  “trapped”  (i.e.,  the 
acceptance  rate  for  particle  transfer  moves  becomes  very  low) 
when  it  acts  as  the  donor  of  a  strong  hydrogen  bond  (as  indicated 
by  a  shortening  of  the  H  O  distance  by  w  0.2  A,  from  around  1.8 
to  about  1.6  A)  to  a  neighboring  water  molecule.  It  should  be 
noted  that  this  trapping  was  only  observed  in  some  BLYP  simula¬ 
tions  and  not  those  performed  with  an  empirical  force  field.  Such 
a  strongly-bound  configuration  may  be  the  precursor  for  a  proton 
transfer  event  but  acid  dissociation  is  not  allowed  in  the  present 
MC  simulations.  In  addition,  different  types  of  local  solvation 
environments  have  also  been  observed  for  undissociated  HN03 
in  aqueous  solution  by  EXAFS  and  first  principles  mole¬ 
cular  dynamics.90  The  EFF  simulations  indicate  that  the  system 
size  effects  in  AGlrans  are  negligible.  The  FPMC  simulations 
exhibit  a  larger  temperature  dependence  with  a  slope  of  about 
20  cal  mob1  K  whereas  the  slope  is  about  one  order  of 
magnitude  smaller  for  the  EFF  simulations. 

Table  2  provides  the  standard-state  Gibbs  free  energy  of  acid 
dissociation  of  HCI  in  aqueous  solution  calculated  with  meth¬ 
ods  SM-ref  and  SM-cluster,  as  well  as  some  of  the  experimental 
data.33  At  T=  300  K  and p  =  1  bar,  both  computational  methods 
predict  AGaaq  that  is  slightly  more  favorable  than  deduced 
from  experiment.  The  SM-ref  and  SM-cluster  methods  yield 
temperature  derivatives  (AGa,aq/A7j  of  10  and  20  cal  mol-1  K-1, 
respectively,  whereas  the  experimental  estimate  is  about 
30  cal  mob1  K  \33  However,  as  discussed  in  the  Introduction, 
there  is  considerable  spread  in  experimental  values  caused  by 
the  difficulty  of  measuring  the  association  constant  of  a  strong 
acid  as  noted  by  Clegg  and  Brimblecombe.24 

The  primary  results  of  this  paper,  namely  AGcomb(7j  for  a 
wide  range  of  temperatures,  are  shown  numerically  in  Table  3 
and  graphically  in  Fig.  1.  For  comparison,  the  experimental 
data  taken  from  Marsh  and  McElroy23  are  also  provided.  We 


Table  2  The  standard-state  Gibbs  free  energy  of  acid  dissociation  (in  kcal  mol-1) 
for  molecular  HCI  and  pKa,aq.The  experimental  values  are  estimated  from  Fig.  1  in 
Pokrovskii33 


T[  K] 

SM-ref 

SM-cluster 

Exp. 

AG°aq 

P^a.aq 

AG°aq 

P/^a,aq 

AG°aq 

300 

-9.7 

-7.08 

-10.2 

-7.45 

-8.4 

350 

-9.3 

-5.83 

-9.2 

-5.75 

-7.2 

400 

-8.8 

-4.83 

-8.2 

-4.47 

-5.5 

450 

-8.3 

-4.04 

-7.1 

-3.46 

-3.8 

Table  3  The  Gibbs  free  energies  of  the  combined  solvation-acid  dissociation 
process  for  HCI  in  water,  in  kcal  mol-1,  calculated  using  eqn  (6)  and  the 
corresponding  Henry's  law  constants.  The  Henry's  law  constants  are  in  units  of 
mol2  kg-2  atm-1 


T[ K] 

Exp.23 

SM-ref 

SM-cluster 

EFF 

BLYP 

EFF 

BLYP 

300 

AUcomb 

-8.6 

-8.9 

N/A 

-9.4 

N/A 

khi  106 

1.86 

3.0 

N/A 

7.3 

N/A 

350 

^^comb 

-7.6 

-8.3 

-7.4 

-8.2 

-7.3 

-Kh/104 

5.97 

16 

4.4 

13 

3.7 

400 

^^comb 

N/A 

-7.7 

-6.0 

-7.0 

-5.3 

Kh/102 

N/A 

170 

18 

71 

7.9 

450 

AUcomb 

N/A 

-7.1 

-4.7 

-6.0 

-3.5 

KhI  10 

N/A 

300 

20 

80 

5.2 

-Z 
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Fig.  1  The  Gibbs  free  energy 

of  the  combined  solvation-acid  dissociation  of  HCI  in 

water  computed  by  eqn  (6),  taken  from  Table  3.  The  experimental,23  BLYP/SM-ref, 
BLYP/SM-cluster,  EFF/SM-ref,  and  EFF/SM-cluster  results  are  depicted  with  X's,  filled 

triangles,  open  triangles,  filled  squares,  and  open  squares, 

respectively. 

feel  justified  in  using  this  data  set  for  comparison  because  it  is 
the  most  complete  data  set  available,  and  questions  were  raised 
primarily  for  the  data  below  temperatures  in  which  we  are 
interested  in  here.  Unfortunately,  no  data  could  be  found  for 
temperatures  exceeding  the  normal  boiling  point  of  water.  The 
agreement  between  experimental  and  computational  predic¬ 
tions  is  very  satisfactory  with  the  predicted  values  at  the  two 
lowest  temperatures  falling  within  1  kcal  mob1  of  the  best 
experimental  estimate.  The  predictions  using  the  BLYP  func¬ 
tional  for  AGtrans  and  either  method  for  AGa  aq  are  even  closer 
to  the  experimental  value  at  T=  350  K.  This  gives  us  confidence 
in  our  predictions  at  higher  temperatures  for  which  experi¬ 
mental  data  are  not  available.  Assuming  that  the  temperature 
dependence  of  the  experimental  data  at  the  lower  temperatures 
is  a  good  guide,  the  BLYP/SM-ref  combination  appears  to  give 
the  most  accurate  predictions.  The  EFF/SM-cluster  combination 
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also  yields  a  good  temperature  dependence,  but  the  AGcomb 
values  are  too  favorable  by  about  0.7  kcal  mol-1.  The  BLYP/ 
SM-cluster  and  EFF/SM-ref  combinations  either  over-  or 
underestimated  the  temperature  dependence  but  agree  quite 
well  with  the  experimental  data  at  ambient  conditions  (using  a 
linear  extrapolation  for  the  BLYP/SM-cluster  method). 

It  can  be  seen  in  Tables  2  and  3  that  the  agreement  between 
the  simulated  and  experimental  results  increases  for  the  combined 
Gibbs  free  energy  compared  to  that  seen  in  the  standard-state 
Gibbs  free  energy  of  acid  dissociation.  This  implies  a  fortuitous 
cancelation  of  errors,  although  the  improvement  is  still  within 
the  estimated  uncertainty  for  each  method.  In  addition,  the 
spread  in  the  calculated  results  increases  with  increasing 
temperatures  (by  about  3  kcal  mol-1  over  100  K). 

For  completeness,  the  KH  values  are  also  listed  in  Table  3 
because  these  allow  for  easier  comparison  with  previous  experi¬ 
mental  data.  As  noted  in  the  Introduction,  Marsh  and  McElroy23 
and  Clegg  and  Brimblecombe24,26’34  report  KH  values  around 
2  x  106  mol2  kg-2  atm-1  at  ambient  conditions.  The  predic¬ 
tions  based  on  EFF  values  for  AGtrans  yield  a  slight  over 
estimation  of  KH  at  T  =  300  K,  and  the  four  computational 
methods  bracket  the  experimental  value  at  350  K. 

4  Summary  and  conclusions 

We  report  the  Gibbs  free  energy  of  the  combined  solvation-acid 
dissociation  of  hydrogen  chloride  in  aqueous  solution;  the  values 
were  obtained  by  using  a  combined  molecular  simulation  and 
implicit  solvation  model  approach.  This  was  done  by  dividing  the 
total  Gibbs  free  energy  change  into  two  distinct  steps.  The  Gibbs  free 
energy  change  in  the  first  step,  transfer  of  molecular  HC1  from  the 
vapor  phase  into  liquid  water,  was  calculated  via  Monte  Carlo 
simulation  in  die  Gibbs  ensemble  using  either  Kohn-Sham  density 
functional  theory  or  an  empirical  force  field  to  describe  the  inter¬ 
actions.  The  second  step,  dissociation  of  molecular  HCl  in  liquid 
water,  was  examined  via  a  continuum  solvation  model  utilizing 
either  experimental  reference  data  or  microsolvated  clusters.  The 
calculations  were  performed  at  four  different  temperatures  ranging 
from  300  to  450  K.  All  four  combined  methods  yield  predictions  with 
satisfactory  accuracy,  but  the  combination  using  the  BLYP  func¬ 
tional  for  the  molecular  simulations  and  the  reference  data  for  the 
solvation  model  stands  out  in  accuracy  with  a  value  of  AGso1v  = 
—7.4  kcal  mol-1  at  T=  350  K  (only  0.2  kcal  mol-1  deviation  from  the 
best  estimate  based  on  various  experimental  data)  and  a  tempera¬ 
ture  dependence  of  about  30  cal  mol-1 1<  '.  With  this  validation  in 
hand,  the  combined  molecular  simulation/implicit  solvation 
approach  is  used  to  predict  experimentally  inaccessible  quantities 
relating  to  the  solvation  of  strong  acids  in  water  and  to  provide 
microscopic-level  information  on  local  solvation  structure.  Such 
solvation  processes  have  wide  applicability  in  chemical  synthesis, 
environmental  chemistry,  and  atmospheric  science. 
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